Tetramerization in a SU(4)-Heisenberg model on the honeycomb lattice 
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The SU(4) Heisenberg model can serve as a low energy model of the Mott insulating state in materials where 
the spins and orbitals are highly symmetric, or in systems of alkaline-earth atoms on optical lattice. Recently, 
it has been argued that on the honeycomb lattice the model exhibits a unique spin-orbital liquid phase with an 
algebraic decay of correlations [P. Corboz et al, Phys. Rev. X 2, 041013 (2012)]. Here we study the instability 
of the algebraic spin-orbital liquid toward spontaneous formation of SU(4) singlet plaquettes (tetramerization). 
Using a variational Monte Carlo approach to evaluate the projected wave-function of fermions with n-fiux state, 
we find that the algebraic liquid is robust, and that a finite value of the next nearest exchange is needed to induce 
tetramerization. We also studied the phase diagram of a model which interpolates between the nearest neighbor 
Heisenberg model and a Hamiltonian for which the singlet-plaquette product state is an exact ground state. 
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PACS numbers: 67.85.-d, 71.10.Fd, 75.10.Jm,75.10.Kt 



I. INTRODUCTION 

In transition metal oxides the partially filled (i-orbitals usu- 
ally undergo a Jahn-Teller distortion which removes the or- 
bital degeneracy. However, in certain cases the orbitals re- 
main fluctuating down to low temperatures at energies com- 
parable to those of spin fluctuations. Possible examples in- 
clude the triangular system LiNiOi [1, 2], and more recently 
Ba3CuSb209 [3, 4] with Cu ions forming a decorated honey- 
comb lattice. Theoretically, the scenario of spin-orbital liquids 
was studied in Refs. [5-12]. 

The minimal model which includes the interaction of spins 
and orbitals is the Kugel-Khomskii model. Due to the spatial 
extension of the orbitals, the Kugel-Khomskii model is usu- 
ally anisotropic in the orbital part. However, under specific 
circumstances, the model acquires its highest symmetry form 
and can be written as 



<K = g(2S,-S, + ij(2T,-T, + ij , 



(1) 



where S,- are spin-1/2 operators acting on the |t) and |X) states 
and T j are operators acting on the twofold degenerate orbital 
degrees of freedom (which we denote by a and b). For this 
special case the system is fully SU(4) symmetric. Denoting 
the local spin-orbital states by colors, !•) =|t a), !•) =|X a), 
!•) -lib), and | ) =\ib), the Hamiltonian (1) is equivalent to 



(2) 
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where P, j, acting in the color space, interchanges the states on 
sites i and j. This Hamiltonian defines the SU(4) symmetric 
Heisenberg model. 

The same model also arises in the Mott insulating state of 
ultra-cold alkaline-earth atoms trapped in an optical lattice. 
In this case, it is the nuclear spin of length F which becomes 
the only relevant degree of freedom with - 2F + 1 states, 
and the interaction leads to the SU(N) symmetric Heisenberg 
model[13-16]. 



The absence of spin and orbital ordering in BasCuSbaOg 
[3, 4] prompted authors of Ref. [17] (including the present 
authors) to study the symmetric Kugel-Khomskii model on 
the honeycomb lattice. Combining the results of a variety of 
numerical and analytical methods, they proposed that an al- 
gebraic spin-orbital liquid (or color liquid) state is realized. 
In this color-liquid the color-color correlation decays alge- 
braically with the distance, and the ground state does not 
break any space group symmetry, nor the SU(4) symmetry. In 
particular, it can be described by a Gutzwiller projected vari- 
ational wave function based on the ;7r-flux state of fermions 
on the honeycomb lattice at 1/4 filling. In the ;r-flux state a 
gauge field corresponding to a magnetic flux of half the flux 
quantum is introduced. The Fermi surface is actually a point, 
as the Fermi-level is at the Dirac point introduced by the n- 
flux. Naturally, the question arises if the variational state is 
stable against the opening of a gap at the Fermi level of the 
unprojected fermions. This can be realized by modulating the 
hoppings and the on-site energies. Here, our aim is to study 
the stability of the algebraic color liquid state against such an 
instability, the tetramerization, using the variational approach. 

The motivation behind choosing tetramerization is not only 
because it opens a gap in the unprojected fermion spectrum, 
but also by the natural tendency of the SU(4) colors to form 
4-site singlet plaquettes [7]. Such plaquettes were observed 
e.g. in two-leg SU(4) ladders [18], and even a Hamiltonian 
based on projection operators can be constructed where de- 
coupled SU(4) singlet plaquettes are the exact ground state 
[19]. Tetramerization also happens in two-dimensional mod- 
els: on the checkerboard lattice the plaquette-product wave 
function is an exact ground state of a suitable Hamiltonian 
which is a sum of projectors, and the singlet formation per- 
sist also in the simple Heisenberg model [20]. For the square 
lattice, the situation is less evident: plaquette formation was 
proposed in Refs. [21-23], but debated in Ref. [10]. Our cur- 
rent understanding is that a dimerization occurs [24]. 

The paper is structured as follows: In Sec. II we construct a 
Hamiltonian which has singlet-product states as exact ground 
states. In Sec. Ill we construct an order parameter to detect 
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tetramerization. In Sec. IV we present results of the vari- 
ational Monte Carlo calculations. Finally we conclude in 
Sec. V. 



II. A MODEL WITH AN EXACT, TETRAMERIZED 
GROUND STATE 

Using projection operators, we can construct a Hamiltonian 
for which a singlet-product wave function is an exact ground 
state, where the four spins of the SU(4) singlets are located 
on 4-site stars (a central site with the 3 nearest neigbor sites). 
There are four such equivalent coverings, that are related by 
translations and inversion. For that reason, we consider the 



(3) 



operator, where P,y is the SU(4) symmetric Heisenberg ex- 
change. For disjunct (ij) and (kl) the = Qlj),(ki) ^o\ds, 
so the Q(ijnki) is a projection operator with eigenvalues and 
1 . It gives if the wave function on the (ij) or the (kl) bonds 
is antisymmetric, i.e. belongs to a singlet. If we choose the 
(ij) and (kl) bonds to be the nearest neighbor parallel bonds, 
as shown in Fig. 1(a), at least one of the (ij) and (kl) bonds is 
covered by a SU(4) singlet, so that 2(/7),«7)l^sp) = for any 
of the four |*I'sp)'s. It follows, that the |^sp) are ground states 
of the 



(/./),(*') 



(4) 



Hamiltonian with energy: since it is a sum of positive def- 
inite operators, the energy of the ground state must be > 0. 
The sum in the 1-{q is over the 3A^ nearest neighbor parallel 
bond configuration. We note that this construction is slightly 
different from those mentioned in Refs. [19, 25, 26]. 




FIG. 1. (Color online) (a) The modulation of the bond energies in 
the tetramerized state on the honeycomb lattice. In the ground state 
of the Hq the spins on the four sites connected by the purple bonds 
make a SU(4) singlet (here emphasized by the magenta triangles), 
these singlet tetramers are decoupled from each other. The two thick 
black bonds represent the configuration of the (i f) and (kl) bonds of 
the Q[ij)xki) projector in the 'Hq Hamiltonian. (b) A SU(4) spin singlet 
that is an eigenstate of the Q projection operators with eigenvalue, 
(c) An invalid SU(4) singlet, as the 2(i5)(78) cannot be satisfied. 

We can also easily convince ourselves that no other type 
of SU(4) singlet coverings satisfies 'Hq. Here we follow the 




FIG. 2. (Color online) Graphical representation of the Ry- = Pi g + 
^2,5 - ^3,8 - Pa,i component of the three-dimensional order param- 
eter: (a) on the lattice (the numbers enumerate equivalent sites in 
the 8-site unit cell), and (b) in the cube. The cube is topologically 
equivalent to ordering in the 8-site unit cell. The bonds partici- 
pating on R^. are shown as double lines. We have assumed that 
one of out of the four tetramerized ground state is realized, in this 
case (Ry,Ry,R-) oc (1, 1, 1). (c) The non-equivalent nearest neighbor 
and next nearest neighbor bond energies in the tetramerized state. 
In the singlet-product wave functions the {P„) = (P,.) = -1 and 
(Pi,) = {Pj) = {Pg) = 1/4, while in the absence of tetramerization 
{Pa) = (Pt) = (P^) and (P,) = (Pj) = {P,) = {P2). 



site numbering shown in Fig. 1: (i) To make Q(i2)(34) satis- 
fied, we choose to antisymmetrize the spins on bond (12); 
(ii) To satisfy 2(i5)(23), we antisymmetrize (23); (iii) To make 
the 2(i7)(26) happy, we need to choose to antysimmetrize spins 
on band (26) [as shown in Fig. 1(b)]. Antisymmetrizing (17) 
would mean that we are not able to satisfy 2(i5)(78) Fig. 1(c), 
as we cannot antisymmetrize colors on more than four sites. 
So, the star shaped singlet in Fig. 1(b) is the only possible 
SU(4) singlet that satisfies all the ^(//irtO- The two choices we 
made in step (i) and (ii) in selecting the bonds of Q on which 
we antisymmetrized the spin lead to the fourfold degeneracy 
of the ground state. 

The following model interpolates between the Heisenberg 
model and the 'Hq: 



(ij) 



(5) 



{ij),(kl) 



For T] - I, the 'H,j = "Hq, while for rj - we recover the 
nearest neighbor SU(4) Heisenberg model, Eq. (2). To find 
the transition point between these two cases, we will study 
this model using the VMC approach in Sec. IV. 



III. TETRAMERIZATION ORDER PARAMETER 

To construct the order parameter of the tetramerization, we 
start from the 8 site unit cell of the tetramerized state shown 
in Fig. 2(a). The connectivity of the sites in the unit cell is 
topologically equivalent to a cube, so the Oh point group of 
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the cube can be used to obtain the order parameter The ex- 
pectation values of the P, ^ exchange operators on the 12 near- 
est neighbor bonds in the unit cell transform according to the 
Aig, Eg, T2g, Tiu, and T2u irreducible representations of the 



Oh group, and can serve to detect order of some kind, 
order parameter of the tetramerization is given by 



Pl.l - P2,S + ^^3.5 - P4,6 
Pl,& - PlJ - Pifi + P4.5 ) 



The 



(6) 



and transforms as the three dimensional T2g irreducible rep- 
resentation. For the singlet product ground state of the "Hq 
shown in Fig. 2(a) R = (-5/2,-5/2,-5/2), while the for 
the remaining three ground states R = (-5/2,5/2,5/2), R = 
(5/2,-5/2,5/2), and R = (5/2,5/2,-5/2). Denoting the 
components of the order parameter as R = {Ry~, R^r, Rxy) [we 
use the (yz, xz, xy) notation as it also transforms according to 
the irreducible representation in O/,], the Landau free en- 
ergy can be written as 



E = -C2 {rI + Ri + R%) + CiRy,R,,R„ + 



(7) 



The cubic invariant with C3 > selects the 4 directions [R oc 
(-1,-1,-1), (-1,1,1), (1,-1,1), and (1,1,-1)], and makes 
the phase transition into the tetramerized phase a first order 
one. 

In the variational calculation, where we explicitly break the 
four-fold degeneracy by hand, it is sufficient to consider the 
difference of the energies on a nearest neighbor bond belong- 
ing to a tetramer and on a nearest neighbor bond connecting 
two tetramers [see Fig. 2(c)]: 



■ ((Pa) - {Ph}) 



(8) 



For the uniform spin-orbital liquid that does not break the 
space group of the lattice r = 0, while for the fully tetramer- 
ized state (Pa) - -1 and {Pi,} - 1/4, so that r = 1 in 
this case. The two order parameters are related as R = 
2«P«>-<Pfc»(l,l,l) = -(5r/2)(l,l,l). 



IV. VARIATIONAL MONTE CARLO RESULTS 

The variational approach is based on the fermionic repre- 
sentation of the SU(N) algebra. The exchange operator Pij 
can be written as 



(9) 



where /j (/) creates a fermion of color a on site ;, and /„ (/) 
annihilates it. To treat the four fermion term, it is custom- 
ary to introduce a bond mean-field decoupling of the P, j [27], 
leading to a free-fermion Hamiltonian, where the hopping am- 
plitudes are determined self-consistently. The mean-field pro- 
cedure often results in solutions where the product of the hop- 
pings around a plaquette is a complex number oc e'*, as if the 




(b) 




FIG. 3. (Color online) (a) The modulation of the hopping amplitudes 
and of the on-site energies in the ;r-fiux hopping Hamiltonian. The 
dashed hoppings represent negative hopping amplitudes. The 8-site 
hexagonal unit cell is shown in orange, (b) The Brillouin zone of the 
honeycomb lattice (blue, with high symmetry points K and M) and 
of the 8-site unit cell (red, with high symmetry points K2 and Mj) 



fermions were picking up a phase due to a magnetic field with 
flux O threading through the plaquette (in convenient units). 
The finite flux can considerably change the band structure. 
Chiral solutions (with <t) 0, tt) were found for the SU(N) 
model with > 5 on the square lattice [28] and for the SU(6) 
model on the honeycomb lattice [29]. While the mean-field 
solutions give a very useful insight into the possible nature of 
the ground state, it describes free fermions where the on-site 
occupancy is fixed on the average only. 

In a complementary approach one does not search for 
mean-field solutions, but takes the ground state wave func- 
tion of some free-fermion Hamiltonian, and by applying a 
Gutzwiller projection one ensures that the occupation of each 
site is exactly one. The projection is done numerically on a 
finite size cluster, using Monte-Carlo importance sampling, 
where one can efliciently sample the wave function and cal- 
culate the energy and correlations [30, 31]. Regarding two- 
dimensional SU(N) Heisenberg models, the method has been 
applied to the SU(4) model on the square [10] and honeycomb 
lattice [17], and to the SU(3) model on the triangular lattice 
[32] and honeycomb lattice [33]. 

As already mentioned in the introduction, the energy of the 
Gutzwiller projected 7r-flux state is much lower than the en- 
ergy of the 0-flux state for the SU(4) model on the honey- 
comb lattice [17]. So to study tetramerization within the vari- 
ational approach, we modify the free fermion Hamiltonian by 
strengthening the hoppings ftet on the bonds belonging to the 
tetramers and by introducing a negative onsite energy Stet for 
the sites at the center of the tetramers, while maintaining the 
;/r-flux state. The Hamiltonian is given by 



(i.j) 



, (10) 



the tij hoppings and e, on-site energies are chosen as shown 
in Fig. 3(a). The sign of the tij is such that the product of 
the hoppings around the hexagons is negative; the product is 
either -f^ or -f^rfef depending on the hexagon. This realizes 
the TT-flux state, the variational state with the lowest energy for 
the nearest neighbor exchange. 
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FIG. 4. (Color online) The two-fold degenerate bands of the free fermions of the tetramerized ;r-flux state for a few selected values of the 
e,et and t,a along a path connecting high-symmetry points in the reduced Brillouin zone [see Fig. 3(b)]. F is the center of the Brillouin-zone 
(k = 0), M2 is the k = (n/i, 0) (and the symmetry related points), and K2 is the k = (0, 2;r/3 VS). Occupied bands are colored dark purple, 
while empty bands are colored light purple. 



A. Free fermions in the A:-space 

Before discussing the results of the Gutzwiller projection, 
it is instructive to look at the properties of free fermions with 
TT-flux per plaquette. Fourier transforming the fermionic op- 
erators in a conveniently chosen gauge, we can write the hop- 
ping hamiltonian in the A:-space as an 8 x 8 matrix of the form 
(following the site ordering shown in Fig. 2) 



-r 
& 



(11) 






2iky f 

c ftet 




hel e 







t 


-t 


" 'tet 


-t 





t 


-iVik.-ikyf 

- 'tet 


t 


-t 






where the only nonzero element of S matrix is fii j - Stet and 



(12) 

describes the hopping between the nearest neighbor sites. In 
the absence of the on-site potential Stet, the £ = and by 
squaring the eigenvalue problem we get 'F^ 'F\tf)} = i.e. 
the bands are particle-hole symmetric. 

Solving for the eigenvalues of Eq. (11), we find four two- 
fold degenerate bands in the Brillouin zone of the 8-site unit 
cell, shown in Fig. 4. The two-fold degeneracy is the direct 
consequence of the tt-Aux. In the limiting case of ftet - 0, 
when we end up with a hexagon and two isolated sites in the 
unit cell, the degeneracy is present as the alternating f hop- 
pings on the hexagon makes the energy levels double degen- 
erate. In the other limit, when t - Q, the situation is sim- 
pler, as there are two decoupled tetramers in the unit cell, 
trivially giving a two-fold degeneracy. To see the degener- 
acy for arbitrary values of hopping parameters, we start from 
a one-fermion wave function where the amplitude is 1 

at the center of one of the tetramers. We can then generate 
a four-dimensional Hilbert-space by successively acting with 



the 'Hf(k) to the \t//iA} state 

l«AlA> 
l«A2A> 



(0,0,0,0,1,0,0,0) 

1 



— (0,1,1,1,0,0,0,0) 
V3 



(13) 
(14) 



1 



0,e 



2ik, 



iy/3k,-ik. 



where 



3Vi-lrkP 

e-'*'-'^*' -rk'0,o,o,o) , 
1 



;(0,0, 



V27Vi-lrkP 

.2*v _ ^-i^k,-ik, J^k,-iky _ ^2ik,^ 



(15) 

e"'^*^«-*> _ gi^k,-ik. 



(16) 



Further acting with '?//(k) on \^aa) will result a linear com- 
bination of Ii/'ia), VI^ia), |iA3a), and \^aa), so that the Hilbert 
space closes. In this Hilbert space the hopping Hamiltonian is 
reduced to 



"~£tet 




V3(i-lrkP)ftet 





Vly^ftet 











V3(i-lrkP)?tet 








V3? 








V3r 






(18) 

The dispersion Sk of the fermions is given by the characteristic 
polynomial of 'H'^ 

= 4- £,et4 - 34 (f2 + 4) + Setetekf' + 9|rkl'f'4 ■ (19) 

Similarly, we can start with where the amplitude is 1 

at site 1 . This leads to an four-dimensional Hilbert space that 
is orthogonal to the {\4i\a), VI^ia), \4''iA), \4'aa)], and we end up 
with a reduced hopping matrix that is the complex conjugate 
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of TYj, with a spectrum also given by Eq. (19). So we have 
succesfully block diagonalized 'Hj and shown that the spec- 
trum is two-fold degenerate. 

In Fig. 4 we show the band structure for different choices 
of the ftet and Stet- The bands of the uniform tt-Aux state of 
the algebraic color liquid are shown in Fig. 4(a). There are 
three Dirac -points in the band structure, two at F point with 
s - ± V3r and one at the M2 point with s = - the latter one 
is counterpart of the Dirac -point of the hopping hamiltonian 
with uniform hoppings, and remains a Dirac point as long as 
Stet = 0, as seen in Fig. 4(b). For the SU(4) problem, however, 
the lower Dirac -point at k = is more important: it makes the 
Fermi surface of the quarter filled Fermi sea point-like and is 
responsible for the algebraic decay of the color-color correla- 
tion functions. Tetramerization removes the Dirac point and 
a gap opens between the lowest band and the band above it 
at the F point. However, fine tuning the parameters, we can 
recover the Dirac point: for k = the Eq. (19) factorizes into 
(e^ - 3f^)(e^ - StetSk - 'ii^ex) - 0. As we change the £,et and 
ftet, a level crossing happens (i.e. the gap closes) at s = - V3f 
when 

£tet = ^ (ft'et - t^) ■ (20) 

A typical spectrum with a Dirac -point for such a particular 
value of etet is shown in Fig. 4(d). We note that in a finite size 
cluster a level crossing happens in the free-fermion spectrum 
at the r point. 



B. Gutzwiller projection 

In the following, we calculate the expectation value of the 
exchange Hamiltonian in the Gutzwiller projected wave func- 
tion using Monte-Carlo importance sampling as we change 
the values of ftet and Stet in the 7r-flux phase. A short descrip- 
tion of the method we use is given in Ref. [17]. We choose 
a 96 site cluster that has the full symmetry of the hexagonal 
lattice and is compatible with the tetramerization pattern. In 
our Monte-Carlo algorhythm an elementary update was the 
exchange of two randomly chosen fermions at arbitrary sites 
with diff'erent colors. Each run had 10'° elementary updates. 
To avoid any autocorrelation effect, we performed measure- 
ments after every 1000 updates, the distance that corresponds 
to around 5 times the autocorrelation length. We averaged 
over all the bonds when calculating the energy, and the errors 
of the energy/site values are of the order of 10"^, thus much 
smaller than the symbol sizes on the plots. Furthermore, the 
boundary conditions (periodic or antiperiodic) of the hopping 
hamiltonian were chosen so that to avoid the degeneracy of 
the free-fermion ground state. 



1. Stability of the Algebraic color liquid 

First, we consider the stability of the algebraic color liquid 
against tetramerization. For that reason we vary the ftet and 
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FIG. 5. (Color online) (a) The energy of the SU(4) nearest neighbor 
Heisenberg model and (b) the order parameter r in the Gutzwiller 
projected variational wave function as a function of the /tct and Stct 
parameters of the hopping Hamiltonian, calculated on a 96 site clus- 
ter. The ftet and Stet were varied in steps of 0. If and 0.2f, respectively. 
Along the dashed black line the free fermion Fermi surface becomes 
a point, following the condition given by Eq. (20). The energy min- 
ima closely follows that line, and similarly the order parameter van- 
ishes in the vicinity of that line. 

Stet in a systematic way and measure the expectation value of 
the energy and the tetramerization order parameter, [Eq. (8)]. 
The energy per site is given by 

E.'JI^llS^, (21) 

where (Pa) and {Pi,} are the expectation value of the exchange 
on the intra- and inter-tetramer bonds, respectively [Fig. 2(c)]. 
Results for the energy as a function of the f,et and Stet are 
shown in Fig. 5(a), together with the tetramerization param- 
eter. Fig. 5(b). It appears that the energy is minimal in a val- 
ley roughly following the condition (20), that is required for 
the Dirac points to form at the Fermi surface, and at the same 
time the tetramerization order parameter disappears. The ac- 
tual minimum is realized for the f^et — f and £tet — 

0: The 

TT-flux State is stable against tetramerization for the nearest 
neighbor exchange on the Heisenberg lattice. This is explic- 
itly seen in Fig. 6(a), where the energy and order parameter 
values of Fig. 5 are replotted with ftet and Stet being implicit 
variables. Quite surprisingly, the tuning of the parameters ftet 
and etet in the free fermion Hamiltonian may have similar ef- 
fects after the Gutzwiller projection. Namely, even for quite 
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FIG. 7. The tetramerization order parameter r as a function of ;/ in 
Eq. (5). A first-order pliase transition occurs at ij 0.5 between 
tlie algebraic color liquid (r = 0) and the tetramerized phase (r > 
0). In the inset we show the expectation value as a function of the 
order parameter close to the phase transition, at = 0.5. Each points 
corresponds to a different value of the free-fermion parameters f,ei 
and Etcti as calculated by the Variational Monte Carlo on the 96 site 
cluster. 



FIG. 6. (Color online) (a) Energy per site of the nearest neighbor 
model, Eq. (21), as a function of the tetramerization parameter r, 
Eq. (8), for various ^tet and £tet values. The singlet product wave 
function is the maximally tetramerized state with r = I and E = 
-9/16 (the rightmost points in the plot). The plot around the r = 
is enlarged in the inset, (b) The expectation values of the exchange 
operators on the non-equivalent first and second neighbor bonds, and 
of the four site PijPki operator. Statistical error is smaller that the 
symbol size. 



different values of ftet and Stet the E vs. r plot essentially lay 
on a single curve. Some deviation from this behavior is seen 
as ftet — * 0, as exemplified by the ftet/f = 0.2 and 0.4 points for 
not too large values of |r| (for larger values of |r| the points for 
these hoppings merge with other points). This behavior is also 
seen for the exchanges on the different non-equivalent bonds 
[Fig. 6(b)]. 



and that the transition can be captured by the variational wave 
function. 

To explore this possibility, we have looked up minimum 
energy solution among the calculated ftet and Stet values for a 
given J], and plotted the corresponding value of order parame- 
ter r. The energy per site is given by 



^ 3 {Pg} + (Pb) 

E = 

2 2 



+ -?l{l+{PijPu}) , 



(22) 



and the result is shown in Fig. 7. The algebraic color liquid 
is stable up to ?/ a; 0.5, where a first order transition happens 
into the tetramerized state - this agrees with the existence of 
the cubic invariant of the order parameter in Eq. (7). The first 
order nature of the transition is further supported by the energy 
versus order parameter plot, shown in the inset of Fig. 7 for 
T] = 0.5; we can clearly recognize the two minima of the en- 
ergy, one corresponding to the algebraic color liquid (r = 0), 
the other at r x 0.3 to the tetramerized state. And, as dis- 
cussed above, for rj - 1, we recover the fully tetramerized 
ground state with r - I. 



2. Phase transition between the algebraic color liquid and the 
tetramerized state 



3. Tetramerization due to next nearest exchange 



For f = the free-fermion Hamiltonian describes decou- 
pled tetramers with localized fermions. The lowest energy 
level of a decoupled star is singly degenerate, and Gutzwiller 
projecting the four fermions ground state gives an SU(4) 
siglet. Not surprisingly, the Gutzwiller projected wave func- 
tion becomes the tetramer-factorized wave function in this 
case, the exact ground state of the of the Hamiltonian (4). So 
changing the 77 in the Hamiltonian (5), we may expect that 
a phase transition occurs at finite value of the 77 parameter 
between the uniform color liquid and the tetramerized state 



Extending the Hamiltonian (2) with the next nearest neigh- 
bor exchange J2, the number of bonds within a singlet 
tetramer increases [c in Fig. 2(c)], but so does the number of 
bonds between the singlets [d and e in Fig. 2(c)] as well. The 
energy per site is now given by 



^ 3/, {Pg) + (Ph) , , (f ,) + 2{Pd} + (Pe) 

E = — h 3/2 



(23) 



where Ji is the nearest neighbor exchange coupling. The 
answer to the question if the J2 drives the system toward 
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FIG. 8. The tetramerization order parameter as a function of JijJi 
in the model with next nearest exchange Ji- A weakly first-order 
phase transition occurs at J2IJ1 ~ 0.38 between the algebraic color 
liquid and the tetramerized phase. In the inset the energy vs. order 
parameter is shown for the calculated /tet and Etct values close to the 
transition point. 



tetramerization is not evident. In the limiting case of the 
tetramer-factorized wave function, the expectation value of 
the energy is £ = -971/16-372/16 = -0.562571-0.187572- 
This energy is to be compared with the expectation value 
of the energy in the algebraic color liquid, where (Pa) - 
(Ph) ~ -0.596 and (Pc) = (Pd) = (Pe) ~ 0.005, so that 
E Si -0.8957i + 0.01572- Comparing these two energies may 
serve as a rough estimate and we learn that tetramerization 
becomes favorable for 72 > 1 .67i . To get a better estimate of 
the transition, we searched for the minimum energy as we in- 
creased the 72, using the tabulated values of the exchange ex- 
pectation for different values of ftet and Stet- In this approach, 
the transition takes place at 72 ~ 0.387i, as seen in Fig. 8. 
The two minima in the energy-order parameter plot (the inset 
in Fig. 8) are much less pronounced than in the case when the 
tetramerization was induced by the four-spin term Q (see the 
corresponding inset in Fig. 7), and, as a consequence, the first 
order character of the phase transition is weaker 



phase transition, as the 72/7i is increased. 



V. CONCLUSIONS 

We studied the tetramerized ground state of the SU(4) 
Heisenberg model on the honeycomb lattice. First, we showed 
that by adding four-site terms of the form PijPu to the Hamil- 
tonian, one can construct a model where the product of four- 
site SU(4) singlets is an exact ground state. This ground state 
breaks the translational invariance and is four-fold degenerate. 
To characterize the tetramerization, we constructed an order 
parameter 

Next, we built a free-fermion Hamiltonian with a ;7r-flux and 
hopping amplitudes ftet and on-site energies £tet that is com- 
patible with the tetramerization pattern. The Gutzwiller pro- 
jected quarter filled Fermi-sea of this Hamiltonian is capable 
of interpolating between the algebraic color liquid of the pure 
Heisenberg model and the tetramer-product wave function of 
the exact ground state for the Hamiltonian with the four-site 
term. 

Using a variational Monte-Carlo method to evaluate the 
Gutzwiller projected fermionic wave function, we calculated 
the energy of the Heisenberg model with the nearest neighbor 
exchange for different values of the hopping parameters ftet 
and on-site energies Stet that controlled the degree of tetramer- 
izaton. In the band structure of the free fermions a Dirac point 
is formed for a particular combination of the ftet and etet- The 
energy is minimal along an extended valley in the (ftet, £tet) 
plane that closely follows the line of the formation of the Dirac 
point, together with a vanishing tetramerization, showing the 
significance of having a Fermi-point instead of a gapped spec- 
trum. Actually, when we plotted the expectation values of the 
exchanges against each other, they roughly fell on single lines. 

The algebraic color-liquid extends in a finite region around 
the nearest neighbor model, and there is a first order phase 
transition between the tetramerized state and the algebraic 
color-liquid as we decreased the strength of the four-site 
terms. Finally, we looked at the effect of the second neigh- 
bor exchange. Within the Variational Monte Carlo method, 
the algebraic color liquid is the ground state for 72 < 0.387i, 
and the tetramerized state wiU have lower energy 72 > 0.387i. 



In the case of very strong 72/7i (i.e. 7i — > 0), we get two 
decoupled triangular lattices. The SU(4) problem on the trian- 
gular lattice was studied in Ref. [34], where a plaquette state 
of resonating SU(4) singlets is proposed as the ground state, 
with a large (12 site) unit cell. Thus we may expect another 
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